# model_age_both
#
# Model Comparative Statics
#

#### Set-Up ####

rm(list = ls())
Sys.info()[7]
# Set file paths 

DATA <- "../external_links/real"
TEMP <- "../temp"
RESULTS <- "../output"

# Go to code folder
setwd("./")
sink("model_comp_stat_xs.txt")


folder <- paste0(TEMP, "/results_exhibits_xs")
if (file.exists(folder)) {
  cat("The folder already exists")
} else {
  dir.create(folder)
}
folder <- paste0(RESULTS, "/results_exhibits_xs")
if (file.exists(folder)) {
  cat("The folder already exists")
} else {
  dir.create(folder)
}

Sys.info()[7]

# Packages
library(MASS)
library(tmvtnorm)
library(haven)
library(dplyr)
library(tidyr)
library(nloptr)
library(corpcor)
library(stargazer)
library(pracma)
library(foreign)
library(writexl)
library(broom)
library(png)
library(tidyverse)
# new packages:
library(compositions) # for rlnorm.rplus


#### Load Data + Assumptions ####

analysis <- read_dta(paste0(DATA, "/model_sample.dta")) %>%
  mutate(oop_i = ifelse(is.na(oop_i)==TRUE, Inf, oop_i))

analysis <- analysis  %>%
  mutate(age_norm = (age - mean(age)) / sd(age))
print(mean(analysis$age))
print(sd(analysis$age))


pa <- 0.005
pfp <- 0.01
pfn <- 0.01

J <- 1 

summary(analysis$fetus_risk)
print(nrow(analysis))
sum(with(analysis, fetus_risk == 2))
analysis <- subset(analysis, fetus_risk != 2)
print(nrow(analysis))
summary(analysis$fetus_risk)


#### Load Necessary Model Functions ####

# source("../estimate_w_xs/r_functions/production_model_xs_new.R")
source("r_functions/production_model_xs_new.R")

# Function that calculates decisions for each observation i 
# NEW: WITHOUT DRAWING A+C (i.e., taking a+c as given)
# NEW: TAKING REC AS GIVEN 
# NEW: ADD POSSIBLE OOP cost FOR INVASIVE = oop2_i 
# NEW: allow mother's prior to be different from the truth (where truth=p_i=KUB)

model_decisions_counterfactuals <- function(par, data, prior){
  
  psi <- par[6]
  
  # NEW: just pull in data (don't draw a+c)
  final <- data
  
  # NEW: allow mother's prior to be different from the truth
  final <- final %>% 
    mutate(mother_prior = (!!prior))
  
  # calculate decisions according to model for each i 
  # recommendation to do invasive based on KUB result causes prior --> prior^psi everywhere on decision tree 
  # NEW: take rec_i as given (i.e., already in data)
  # NEW: use prior (mother_prior), not truth (p_i)
  final <- final %>% 
    mutate(p = ifelse(rec_i==1, mother_prior^psi, mother_prior))
  
  final <- final %>% 
    mutate(true_pr_pos = (1-p_i)*pfp + p_i*(1-pfn) ) %>% # note: uses truth=p_i, not p (which is based on prior)
    mutate(pr_pos = (1-p)*pfp + p*(1-pfn) ) %>%
    mutate(pr_c_pos = (1-pfn)*p/pr_pos ) %>%
    mutate(pr_c_neg = pfn*p/(1-pr_pos) ) 
  
  final <- final %>% 
    mutate(ui2_Y0 = pa*a_i + (1-pa)*p*pmax(a_i, c_i) - oop2_i ) %>% 
    mutate(ui2_N0 = pmax(a_i, p*c_i) ) %>% 
    mutate(ui2_YP = pa*a_i + (1-pa)*pr_c_pos*pmax(a_i, c_i) - oop_i - oop2_i ) %>% 
    mutate(ui2_NP = pmax(a_i, pr_c_pos*c_i) - oop_i ) %>% 
    mutate(ui2_YN = pa*a_i + (1-pa)*pr_c_neg*pmax(a_i, c_i) - oop_i - oop2_i ) %>% 
    mutate(ui2_NN = pmax(a_i, pr_c_neg*c_i) - oop_i ) 
  
  final <- final %>%
    mutate(invasive_P = ifelse(round(ui2_YP, 5) >= round(ui2_NP, 5), 1, 0)) %>%
    mutate(invasive_N = ifelse(round(ui2_YN, 5) >= round(ui2_NN, 5), 1, 0)) %>%
    mutate(invasive_0 = ifelse(round(ui2_Y0, 5) >= round(ui2_N0, 5), 1, 0)) 
  
  final <- final %>%
    mutate(ui1_Y = pr_pos*pmax(ui2_YP, ui2_NP) + (1-pr_pos)*pmax(ui2_YN, ui2_NN) ) %>%
    mutate(ui1_N = pmax(ui2_Y0, ui2_N0) )
  
  final <- final %>%
    mutate(nipt_indiff = ifelse(round(ui1_Y, 5) == round(ui1_N, 5), 1, 0)) %>%
    mutate(pred_nipt = ifelse(nipt_indiff==0, 
                              ifelse(ui1_Y > ui1_N, 1, 0), 0))
  
  final <- final %>%
    mutate(pred_invasive = ifelse(pred_nipt==0, invasive_0, ifelse(sim_nipt_result == 1, invasive_P, invasive_N))) 
  
  # WTP is Ui1_Y - Ui1_N (without factoring in oop_i), need to write it out to get rid of oop_i from eqs
  final <- final %>% 
    mutate(ui2_Y0_wtp = pa*a_i + (1-pa)*p*pmax(a_i, c_i) - oop2_i ) %>% 
    mutate(ui2_N0_wtp = pmax(a_i, p*c_i) ) %>% 
    mutate(ui2_YP_wtp = pa*a_i + (1-pa)*pr_c_pos*pmax(a_i, c_i) - oop2_i ) %>% 
    mutate(ui2_NP_wtp = pmax(a_i, pr_c_pos*c_i)) %>% 
    mutate(ui2_YN_wtp = pa*a_i + (1-pa)*pr_c_neg*pmax(a_i, c_i) - oop2_i ) %>% 
    mutate(ui2_NN_wtp = pmax(a_i, pr_c_neg*c_i)) 
  final <- final %>% 
    mutate(wtp_nipt = pr_pos*pmax(ui2_YP_wtp, ui2_NP_wtp) + (1-pr_pos)*pmax(ui2_YN_wtp, ui2_NN_wtp) 
           - ui1_N)
  
  final <- final %>% #note: the following is a bit diff
    dplyr::select(pregnancy, did_nipt, did_invasive, live_birth, p_i, wave, oop_i, age, sim_positive, 
                  a_i, c_i, oop2_i, rec_i, ave_kub_age,
                  pred_nipt, pred_invasive, invasive_P, invasive_N, sim_nipt_result, wtp_nipt)  
  
  return(final)
}


#### Expand data 1000x and draw a_i and c_i  ####

# Pull in estimated parameters
f_a = ~ 0 + age_norm
f_c = ~ 0 + age_norm
x_vars <- c("age_norm")

unique_flag <- 1

print("Estimated parameters")

est_par <- read.csv(paste0(RESULTS, "/estimate_w_xs/results_no_xs.csv")) %>%
  as.matrix() 
est_par <- est_par[-9:-10,]
print(est_par)


#calculate true positive indicator for chrom. ab. = live birth with chrom. ab. (using actual, not simulated, data), and simulate positive indicator for counterfactuals
# actual <- analysis %>%
#   mutate(act_positive = ifelse(chroma == 1, 1, 0))
actual <- analysis

#Step 0a: expand data K times
K <- 1000
actual_1 <- actual
if(K > 1){
  for (i in 1:(K-1)) {
    actual <- rbind(actual, actual_1)
  }
}
rm(i, actual_1)
actual <- actual %>% arrange(pregnancy)

#Step 0b: simulate positive indicator for chromAb 
actual <- actual %>% 
  mutate(draw = runif(nrow(.), min=0, max=1)) %>%
  mutate(sim_positive = ifelse(draw <= p_i, 1, 0)) %>% 
  select(-draw)

#Step 0c: simulate result of NIPT test, if taken (using simulated chromAb indicator)
actual <- actual %>% 
  mutate(draw = runif(nrow(.), min=0, max=1)) %>%
  mutate(sim_nipt_result = ifelse(sim_positive == 1, 
                                  ifelse(draw > pfn, 1, 0),
                                  ifelse(draw <= pfp, 1, 0))) %>% 
  select(-draw)

# Step 1: draw (a_i, c_i) (using parameter estimates)
set.seed(101)
simulated <- as.data.frame(draw_a_c(est_par, actual, f_a, f_c, unique_flag, x_vars))
nrow(simulated)
simulated <- simulated %>% dplyr::select(a_i, c_i)
# simulated <- cbind(actual, simulated) %>% 
#   dplyr::select(pregnancy, did_nipt, did_invasive, live_birth, p_i, wave, oop_i, age, ave_kub_age, sim_positive, act_positive, a_i, c_i, sim_nipt_result)
simulated <- cbind(actual, simulated) %>% 
  dplyr::select(pregnancy, did_nipt, did_invasive, live_birth, p_i, wave, oop_i, age, ave_kub_age, sim_positive, a_i, c_i, sim_nipt_result)
rm(actual)

# Step 2: Get values of the percentiles to use for comparative static simulated decisions 
# NEW: 6/16/22 add 75th percentile
a_percentiles<- quantile(simulated$a_i, c(0.10, 0.5, 0.75, 0.9))
c_percentiles<- quantile(simulated$c_i, c(0.10, 0.5, 0.75, 0.9))

# Set up results datasets
comp_model <- simulated %>% dplyr::select(pregnancy, p_i, wave)
comp_51_200 <- simulated %>% dplyr::select(pregnancy, p_i, wave)
comp_hr_cov <- simulated %>% dplyr::select(pregnancy, p_i, wave)


# With drawn a_i and c_i 
comp_data <- simulated %>%
  mutate(oop2_i = 0) %>%
  mutate(rec_i = ifelse(p_i >= 1/200, 1, 0))
dec <- model_decisions_counterfactuals(par = est_par, data = comp_data, prior = quo(p_i)) %>% 
  dplyr::select(pred_nipt, pred_invasive, wtp_nipt)

dec <- dec %>%
  rename(wtp_nipt = wtp_nipt) %>%
  rename(pred_nipt = pred_nipt) %>% 
  rename(pred_inv = pred_invasive) %>%
  mutate(pred_invN = ifelse(pred_nipt == 0, pred_inv, NA)) %>%
  mutate(pred_invY = ifelse(pred_nipt == 1, pred_inv, NA))
comp_model <- cbind(comp_model, dec)


## Comparative statics over a_i ##
count <- 0
for (a_fix in a_percentiles) {
  count <- count + 1
  comp_data <- simulated %>% mutate(a_i = a_fix)
  print(count)
  print(a_fix)
  
  # Under the actual rules
  comp_data <- comp_data %>%
    mutate(oop2_i = 0) %>%
    mutate(rec_i = ifelse(p_i >= 1/200, 1, 0))
  dec <- model_decisions_counterfactuals(par = est_par, data = comp_data, prior = quo(p_i)) %>% 
    dplyr::select(pred_nipt, pred_invasive, wtp_nipt)
  
  if (count == 1) {
    dec <- dec %>% 
      rename(wtp_nipt_a_10 = wtp_nipt) %>%
      rename(pred_nipt_a_10 = pred_nipt) %>% 
      rename(pred_inv_a_10 = pred_invasive) %>%
      mutate(pred_invN_a_10 = ifelse(pred_nipt_a_10 == 0, pred_inv_a_10, NA)) %>%
      mutate(pred_invY_a_10 = ifelse(pred_nipt_a_10 == 1, pred_inv_a_10, NA))
  }
  if (count == 2) {
    dec <- dec %>%       
      rename(wtp_nipt_a_50 = wtp_nipt) %>%
      rename(pred_nipt_a_50 = pred_nipt) %>% 
      rename(pred_inv_a_50 = pred_invasive) %>%
      mutate(pred_invN_a_50 = ifelse(pred_nipt_a_50 == 0, pred_inv_a_50, NA)) %>%
      mutate(pred_invY_a_50 = ifelse(pred_nipt_a_50 == 1, pred_inv_a_50, NA))
  }
  if (count == 3) {
    dec <- dec %>% 
      rename(wtp_nipt_a_75 = wtp_nipt) %>%
      rename(pred_nipt_a_75 = pred_nipt) %>% 
      rename(pred_inv_a_75 = pred_invasive) %>%
      mutate(pred_invN_a_75 = ifelse(pred_nipt_a_75 == 0, pred_inv_a_75, NA)) %>%
      mutate(pred_invY_a_75 = ifelse(pred_nipt_a_75 == 1, pred_inv_a_75, NA))
  }
  if (count == 4) {
    dec <- dec %>% 
      rename(wtp_nipt_a_90 = wtp_nipt) %>%
      rename(pred_nipt_a_90 = pred_nipt) %>% 
      rename(pred_inv_a_90 = pred_invasive) %>%
      mutate(pred_invN_a_90 = ifelse(pred_nipt_a_90 == 0, pred_inv_a_90, NA)) %>%
      mutate(pred_invY_a_90 = ifelse(pred_nipt_a_90 == 1, pred_inv_a_90, NA))
  }
  comp_model <- cbind(comp_model, dec)
  
  
  # Under 1/51 to 1/200 coverage rules
  comp_data <- comp_data %>%
    mutate(oop_i = ifelse((p_i >= 1/200) & (p_i <= 1/51), 0, 567.50)) %>% 
    mutate(oop2_i = 0) %>%
    mutate(rec_i = ifelse(p_i >= 1/200, 1, 0))
  dec <- model_decisions_counterfactuals(par = est_par, data = comp_data, prior = quo(p_i)) %>% 
    dplyr::select(pred_nipt, pred_invasive)
  
  if (count == 1) {
    dec <- dec %>% 
      rename(pred_nipt_a_10 = pred_nipt) %>% 
      rename(pred_inv_a_10 = pred_invasive) %>%
      mutate(pred_invN_a_10 = ifelse(pred_nipt_a_10 == 0, pred_inv_a_10, NA)) %>%
      mutate(pred_invY_a_10 = ifelse(pred_nipt_a_10 == 1, pred_inv_a_10, NA))
  }
  if (count == 2) {
    dec <- dec %>% 
      rename(pred_nipt_a_50 = pred_nipt) %>% 
      rename(pred_inv_a_50 = pred_invasive) %>%
      mutate(pred_invN_a_50 = ifelse(pred_nipt_a_50 == 0, pred_inv_a_50, NA)) %>%
      mutate(pred_invY_a_50 = ifelse(pred_nipt_a_50 == 1, pred_inv_a_50, NA))
  }
  if (count == 3) {
    dec <- dec %>% 
      rename(pred_nipt_a_75 = pred_nipt) %>% 
      rename(pred_inv_a_75 = pred_invasive) %>%
      mutate(pred_invN_a_75 = ifelse(pred_nipt_a_75 == 0, pred_inv_a_75, NA)) %>%
      mutate(pred_invY_a_75 = ifelse(pred_nipt_a_75 == 1, pred_inv_a_75, NA))
  }
  if (count == 4) {
    dec <- dec %>% 
      rename(pred_nipt_a_90 = pred_nipt) %>% 
      rename(pred_inv_a_90 = pred_invasive) %>%
      mutate(pred_invN_a_90 = ifelse(pred_nipt_a_90 == 0, pred_inv_a_90, NA)) %>%
      mutate(pred_invY_a_90 = ifelse(pred_nipt_a_90 == 1, pred_inv_a_90, NA))
  }
  comp_51_200 <- cbind(comp_51_200, dec)
  
  
  # Under 1/2 to 1/200 coverage rules
  comp_data <- comp_data %>%
    mutate(oop_i = ifelse((p_i >= 1/200), 0, 567.50)) %>% 
    mutate(oop2_i = 0) %>%
    mutate(rec_i = ifelse(p_i >= 1/200, 1, 0))
  dec <- model_decisions_counterfactuals(par = est_par, data = comp_data, prior = quo(p_i)) %>% 
    dplyr::select(pred_nipt, pred_invasive)
  if (count == 1) {
    dec <- dec %>% 
      rename(pred_nipt_a_10 = pred_nipt) %>% 
      rename(pred_inv_a_10 = pred_invasive) %>%
      mutate(pred_invN_a_10 = ifelse(pred_nipt_a_10 == 0, pred_inv_a_10, NA)) %>%
      mutate(pred_invY_a_10 = ifelse(pred_nipt_a_10 == 1, pred_inv_a_10, NA))
  }
  if (count == 2) {
    dec <- dec %>% 
      rename(pred_nipt_a_50 = pred_nipt) %>% 
      rename(pred_inv_a_50 = pred_invasive) %>%
      mutate(pred_invN_a_50 = ifelse(pred_nipt_a_50 == 0, pred_inv_a_50, NA)) %>%
      mutate(pred_invY_a_50 = ifelse(pred_nipt_a_50 == 1, pred_inv_a_50, NA))
  }
  if (count == 3) {
    dec <- dec %>% 
      rename(pred_nipt_a_75 = pred_nipt) %>% 
      rename(pred_inv_a_75 = pred_invasive) %>%
      mutate(pred_invN_a_75 = ifelse(pred_nipt_a_75 == 0, pred_inv_a_75, NA)) %>%
      mutate(pred_invY_a_75 = ifelse(pred_nipt_a_75 == 1, pred_inv_a_75, NA))
  }
  if (count == 4) {
    dec <- dec %>% 
      rename(pred_nipt_a_90 = pred_nipt) %>% 
      rename(pred_inv_a_90 = pred_invasive) %>%
      mutate(pred_invN_a_90 = ifelse(pred_nipt_a_90 == 0, pred_inv_a_90, NA)) %>%
      mutate(pred_invY_a_90 = ifelse(pred_nipt_a_90 == 1, pred_inv_a_90, NA))
  }
  comp_hr_cov <- cbind(comp_hr_cov, dec)
}

## Comparative statics over c_i ##
count <- 0
for (c_fix in c_percentiles) {
  count <- count + 1
  comp_data <- simulated %>% mutate(c_i = c_fix)
  print(count)
  print(c_fix)
  
  
  # Under the actual rules
  comp_data <- comp_data %>%
    mutate(oop2_i = 0) %>%
    mutate(rec_i = ifelse(p_i >= 1/200, 1, 0))
  dec <- model_decisions_counterfactuals(par = est_par, data = comp_data, prior = quo(p_i)) %>% 
    dplyr::select(pred_nipt, pred_invasive, wtp_nipt)
  if (count == 1) {
    dec <- dec %>% 
      rename(wtp_nipt_c_10 = wtp_nipt) %>%
      rename(pred_nipt_c_10 = pred_nipt) %>% 
      rename(pred_inv_c_10 = pred_invasive) %>%
      mutate(pred_invN_c_10 = ifelse(pred_nipt_c_10 == 0, pred_inv_c_10, NA)) %>%
      mutate(pred_invY_c_10 = ifelse(pred_nipt_c_10 == 1, pred_inv_c_10, NA))
  }
  if (count == 2) {
    dec <- dec %>% 
      rename(wtp_nipt_c_50 = wtp_nipt) %>%
      rename(pred_nipt_c_50 = pred_nipt) %>% 
      rename(pred_inv_c_50 = pred_invasive) %>%
      mutate(pred_invN_c_50 = ifelse(pred_nipt_c_50 == 0, pred_inv_c_50, NA)) %>%
      mutate(pred_invY_c_50 = ifelse(pred_nipt_c_50 == 1, pred_inv_c_50, NA))
  }
  if (count == 3) {
    dec <- dec %>% 
      rename(wtp_nipt_c_75 = wtp_nipt) %>%
      rename(pred_nipt_c_75 = pred_nipt) %>% 
      rename(pred_inv_c_75 = pred_invasive) %>%
      mutate(pred_invN_c_75 = ifelse(pred_nipt_c_75 == 0, pred_inv_c_75, NA)) %>%
      mutate(pred_invY_c_75 = ifelse(pred_nipt_c_75 == 1, pred_inv_c_75, NA))
  }
  if (count == 4) {
    dec <- dec %>% 
      rename(wtp_nipt_c_90 = wtp_nipt) %>%
      rename(pred_nipt_c_90 = pred_nipt) %>% 
      rename(pred_inv_c_90 = pred_invasive) %>%
      mutate(pred_invN_c_90 = ifelse(pred_nipt_c_90 == 0, pred_inv_c_90, NA)) %>%
      mutate(pred_invY_c_90 = ifelse(pred_nipt_c_90 == 1, pred_inv_c_90, NA))
  }
  comp_model <- cbind(comp_model, dec)
  
  # Under 1/51 to 1/200 coverage rules
  comp_data <- comp_data %>%
    mutate(oop_i = ifelse((p_i >= 1/200) & (p_i <= 1/51), 0, 567.50)) %>% 
    mutate(oop2_i = 0) %>%
    mutate(rec_i = ifelse(p_i >= 1/200, 1, 0))
  dec <- model_decisions_counterfactuals(par = est_par, data = comp_data, prior = quo(p_i)) %>% 
    dplyr::select(pred_nipt, pred_invasive)
  if (count == 1) {
    dec <- dec %>% 
      rename(pred_nipt_c_10 = pred_nipt) %>% 
      rename(pred_inv_c_10 = pred_invasive) %>%
      mutate(pred_invN_c_10 = ifelse(pred_nipt_c_10 == 0, pred_inv_c_10, NA)) %>%
      mutate(pred_invY_c_10 = ifelse(pred_nipt_c_10 == 1, pred_inv_c_10, NA))
  }
  if (count == 2) {
    dec <- dec %>% 
      rename(pred_nipt_c_50 = pred_nipt) %>% 
      rename(pred_inv_c_50 = pred_invasive) %>%
      mutate(pred_invN_c_50 = ifelse(pred_nipt_c_50 == 0, pred_inv_c_50, NA)) %>%
      mutate(pred_invY_c_50 = ifelse(pred_nipt_c_50 == 1, pred_inv_c_50, NA))
  }
  if (count == 3) {
    dec <- dec %>% 
      rename(pred_nipt_c_75 = pred_nipt) %>% 
      rename(pred_inv_c_75 = pred_invasive) %>%
      mutate(pred_invN_c_75 = ifelse(pred_nipt_c_75 == 0, pred_inv_c_75, NA)) %>%
      mutate(pred_invY_c_75 = ifelse(pred_nipt_c_75 == 1, pred_inv_c_75, NA))
  }
  if (count == 4) {
    dec <- dec %>% 
      rename(pred_nipt_c_90 = pred_nipt) %>% 
      rename(pred_inv_c_90 = pred_invasive) %>%
      mutate(pred_invN_c_90 = ifelse(pred_nipt_c_90 == 0, pred_inv_c_90, NA)) %>%
      mutate(pred_invY_c_90 = ifelse(pred_nipt_c_90 == 1, pred_inv_c_90, NA))
  }
  comp_51_200 <- cbind(comp_51_200, dec)
  
  
  # Under 1/2 to 1/200 coverage rules
  comp_data <- comp_data %>%
    mutate(oop_i = ifelse((p_i >= 1/200), 0, 567.50)) %>% 
    mutate(oop2_i = 0) %>%
    mutate(rec_i = ifelse(p_i >= 1/200, 1, 0))
  dec <- model_decisions_counterfactuals(par = est_par, data = comp_data, prior = quo(p_i)) %>% 
    dplyr::select(pred_nipt, pred_invasive)
  if (count == 1) {
    dec <- dec %>% 
      rename(pred_nipt_c_10 = pred_nipt) %>% 
      rename(pred_inv_c_10 = pred_invasive) %>%
      mutate(pred_invN_c_10 = ifelse(pred_nipt_c_10 == 0, pred_inv_c_10, NA)) %>%
      mutate(pred_invY_c_10 = ifelse(pred_nipt_c_10 == 1, pred_inv_c_10, NA))
  }
  if (count == 2) {
    dec <- dec %>% 
      rename(pred_nipt_c_50 = pred_nipt) %>% 
      rename(pred_inv_c_50 = pred_invasive) %>%
      mutate(pred_invN_c_50 = ifelse(pred_nipt_c_50 == 0, pred_inv_c_50, NA)) %>%
      mutate(pred_invY_c_50 = ifelse(pred_nipt_c_50 == 1, pred_inv_c_50, NA))
  }
  if (count == 3) {
    dec <- dec %>% 
      rename(pred_nipt_c_75 = pred_nipt) %>% 
      rename(pred_inv_c_75 = pred_invasive) %>%
      mutate(pred_invN_c_75 = ifelse(pred_nipt_c_75 == 0, pred_inv_c_75, NA)) %>%
      mutate(pred_invY_c_75 = ifelse(pred_nipt_c_75 == 1, pred_inv_c_75, NA))
  }
  if (count == 4) {
    dec <- dec %>% 
      rename(pred_nipt_c_90 = pred_nipt) %>% 
      rename(pred_inv_c_90 = pred_invasive) %>%
      mutate(pred_invN_c_90 = ifelse(pred_nipt_c_90 == 0, pred_inv_c_90, NA)) %>%
      mutate(pred_invY_c_90 = ifelse(pred_nipt_c_90 == 1, pred_inv_c_90, NA))
  }
  comp_hr_cov <- cbind(comp_hr_cov, dec)
}

# Save datasets with comp_statics to stata format for plotting
write.dta(comp_model, paste0(TEMP, "/results_exhibits_xs/comp_model.dta"))

## HOW THIS IS CODED: For each ieration within set of coverage rules, save the relevant testing variables (pred_nipt, pred_invasive, pred_invN, pred_invY) and 
# name them eg pred_nipt_a_10 for the iteration with a at 10th percentile. Then, within each coverage regime, create big dataset with p_i and all the testing vars
# across the different comparative statics. Then export this dataset to stata and do the plotting in stata. 




